#include "CompBioTools.hh"
#include "MerHashMap.hh"


RandomNumberGenerator RNGX;

double gIntRange;

void GetRandomMer(Mer &M)
{
  M.mUpper32 = (unsigned int) RNGX.Uniform(0,gIntRange);
  M.mLower32 = (unsigned int) RNGX.Uniform(0,gIntRange);
}


void main(int argc,char** argv)
{

  int subset = atoi(argv[2]);

  izipstream fin;
  fin.open(argv[1]);

  Mer2IntMap SourceMap;
  if (fin.fail()){
	 cerr<<"Couldn't open "<<argv[1]<<endl;
  }
  cout<<"Adding "<<argv[1]<<" to source hash..."<<flush;
  if (subset == -1){
	 SourceMap.Add(fin);
  }else{
	 SourceMap.Add(fin,subset);
  }
  cout<<"done."<<endl;

  cout<<"Source size="<<SourceMap.size()<<" TotalMers="<<SourceMap.mTotalMersInHash<<endl;

  Mer2IntMap::iterator it;
  //for(it = SourceMap.begin();it != SourceMap.end();it++){
  //		Mer temp = (*it).first;
  //		Sequence S;
  //		temp.ToSequence(S);
  //		cout<<S<<" has count "<<(*it).second<<endl;
  // }

  MerSet Mers;
  Mer2IntMap MerMap;
 
  int idx = 0;
  for(it = SourceMap.begin();it != SourceMap.end();it++){

	 Mer M = (*it).first;
	 
	 idx ++;

	 MerMap.Add(M);
	 Mers.Add(M);   // Keep track of all mers, fwd or reverse. 

	 Sequence Seq;
	 M.ToSequence(Seq);

	 if (!(idx%10000)) {
		cout<<"Seq = "<<Seq<<"\t";
	 }

	 Seq.ReverseComplement();
	 if (!(idx%10000)) {
		cout<<"REVSeq = "<<Seq<<endl;
	 }

	 Mer RCM(Seq);

	 MerMap.Add(RCM);
	 Mers.Add(RCM);   // Keep track of all mers, fwd or reverse. 
  }

  // Add a few again so that we will have some count=3 to look at. 
  idx = 0;
  for(it = SourceMap.begin();it != SourceMap.end();it++){
	 
	 Mer M = (*it).first;
	 MerMap.Add(M);

	 idx++;
	 if (idx > 50) break;
  }

  // cout<<"Record RAM"<<endl;
  //string ans;
  //cin>>ans;

  cout<<"Checking that all mers are in Mer2IntMap..."<<endl;
  MerSet::iterator sit;
  for(sit = Mers.begin();sit != Mers.end();sit++){
	 Mer testMer = (*sit);
	 if (!MerMap.Contains(testMer)){
		Sequence S;
		 testMer.ToSequence(S);
		cout<<S<<"\t is not in Mer2IntMap"<<endl;
	 }
  }
  cout<<"done."<<endl;

  cout<<"Checking that all mers have count =2 ..."<<flush;
  for(it = MerMap.begin();it!= MerMap.end();it++){
	 Mer temp = (*it).first;	 
	 Sequence S;
	 temp.ToSequence(S);
	 if (!SourceMap.Contains(temp)){
		cout<<"Source Map does not contain "<<S<<endl;
	 }else{
		int sourceCount = SourceMap[temp];
		//if (sourceCount != (*it).second){
		if ((*it).second != 2){
		  cout<<S<<" has count "<<(*it).second<<" where source has count "<<sourceCount<<endl;
		}
		  //}
	 }	
  }
  cout<<"done."<<endl;

  cout<<"Mers.size()="<<Mers.size()<<endl;
 

  cout<<"MerMap.size()="<<MerMap.size()<<endl;
  cout<<"MerMap.total="<<MerMap.mTotalMersInHash<<endl;
  
  
}


// Entropy tests with symmetric Mer2IntMap
// 
//   2,000,000 input mers  25.1M   12.55 bytes/input mer  25 bytes/symm mer
//   8,000,000 input mers  91.4M   11.42 bytes/input mer  23 bytes/symm mer
// 
// So, based on this, should be able to get about 1*10^9/(23) = 43 million symm mers
// or 86 million raw mers.   There are 6*10^9 possible raw mers, so 69.7 batches
// should be needed.  
// 
// Batch based on 7 random bits of mUpper32, so there should be 128 batches.  
// Should take about 4 cluster hours.  



// Woodstock tests. 
//
//   2,000,000 input mers (fwd+rev)  22M   11 bytes/mer   
//  20,000,000 input mers (fwd+rev) 201M   10 bytes/mer 
//  80,000,000 input mers (fwd+rev) 803M   10 bytes/mer 
// 100,000,000 input mers (fwd+rev) 
